/*
 Class SetPart
 Description: exact set partition model modified by Pablo Munhoz
 Author: Sabir Ribas
 Email: sabirribas@gmail.com
 Date: Jun 22, 2010
 Version: 0.6
 */

#ifndef SETPART_HPP
#define SETPART_HPP

#include <cstdio>
#include <cstdlib>
#include <iostream>
#include <sstream>
#include <vector>

#include <ilcplex/ilocplex.h>
//ILOSTLBEGIN

#include "../Heuristics/BasicIteratedLocalSearchPRVP.hpp"

struct cutinfo
{
	CPXLPptr model;
	int numcols;
	int N;
	int P;//number of periods
	double Q;
	double *delivery;
	double *x;
	int numeroDeRotas;
	vector<vector<int> > rotasDoCliente;
	vector<vector<int> > rotas;
	vector<int> periodos;
	vector<vector<int> >* dados;
};

typedef struct cutinfo CUTINFO, *CUTINFOptr;

struct incumbentinfo
{
	CPXLPptr model;
	int numcols;
	BasicIteratedLocalSearchPRVP* ils;
	int P; //numero de periodos
	vector<vector<int> > rotasDoCliente;
	vector<vector<int> > rotas;
	vector<int> periodos;
	PRVPProblemInstance* pPRVP;
	double tolerancia;
};

typedef struct incumbentinfo INCUMBENTINFO, *INCUMBENTINFOptr;

double custoBestOfAll;
double cutOff;
SolutionPRVP* solution;
EvaluationPRVP* evaluation;

#define EpsForIntegrality 0.000001

using namespace std;

class SetPart
{
protected:
	PRVPProblemInstance& pPRVP;

	// Funcao atualizaCutinfo

	void atualizaCutinfo(CUTINFO* cutinfo, CPXLPptr model, LPool& lista)
	{
		int nCli = pPRVP.nCli;
		cutinfo->model = model;
		cutinfo->numcols = lista.rotas.size();
		cutinfo->N = nCli;
		cutinfo->P = pPRVP.nPeriod;
		cutinfo->Q = pPRVP.capacidade;
		cutinfo->x = new double[lista.rotas.size()];
		cutinfo->delivery = new double[nCli + 1];
		cutinfo->numeroDeRotas = lista.rotas.size();
		cutinfo->rotasDoCliente = *lista.cliInRoute;
		cutinfo->rotas = lista.rotas;
		cutinfo->periodos = lista.periodo;
		cutinfo->dados = pPRVP.dados;

		for (int i = 0; i < nCli + 1; i++)
		{
			cutinfo->delivery[i] = pPRVP.dados->at(i).at(DEMANDA);
		}
	}

	void atualizaIncumbentinfo(INCUMBENTINFO* incumbentinfo, CPXLPptr model, BasicIteratedLocalSearchPRVP& ils, LPool& lista)
	{
		incumbentinfo->model = model;
		incumbentinfo->numcols = lista.rotas.size();
		incumbentinfo->ils = &ils;
		incumbentinfo->rotasDoCliente = *lista.cliInRoute;
		incumbentinfo->rotas = lista.rotas;
		incumbentinfo->periodos = lista.periodo;
		incumbentinfo->P = pPRVP.nPeriod;
		incumbentinfo->pPRVP = &pPRVP;
		incumbentinfo->tolerancia = lista.tolerancia;
	}

	static string getCollumName(CPXCENVptr cpxEnv, CPXLPptr cpxModel, int collumIndex)
	{
		int surplus;
		CPXgetcolname(cpxEnv, cpxModel, NULL, NULL, 0, &surplus, collumIndex, collumIndex);
		int storespace = -surplus;

		char* aux = new char[storespace];
		char** colname;
		colname = new char*[1];
		//colname[0] = new char[largestNameLength];
		strcat(colname[0], "");

		int largestNameLength = 100;
		CPXgetcolname(cpxEnv, cpxModel, colname, aux, largestNameLength, &surplus, collumIndex, collumIndex);

		string result(colname[0]);

		//delete[] colname[0];
		delete[] colname;
		delete[] aux;
		return result;
	}

	static void storeLPSolution(CPXCENVptr env, CPXLPptr model, int numcols, double *x, double *sol)
	{
		char name[100];
		for (int k = 0; k < numcols; k++)
		{
			string s = SetPart::getCollumName(env, model, k);
			int i;
			if (s.at(0) == 'Y')
			{

				strncpy(name, s.c_str(), 100);
				sscanf(name, "Y(%d)", &i);
				sol[i] = x[k];
			}
		}
	}

	static bool estaPadrao(CUTINFOptr cutinfo, vector<int>& encontrados, int c)
	{
		bool correto = false;
		int numComb = cutinfo->dados->at(c).at(NCOMB);
		//cout << "client " << c << endl;
		//cout << "encontrado: " << encontrados << endl;
		for (int i = 0; i < numComb; i++)
		{
			const vector<int>& escolha = PRVPProblemInstance::decompoeEscolha(cutinfo->dados->at(c).at(NCOMB + 1 + i), cutinfo->P);
			//cout << "escolha [" << i << "]" << escolha << endl;
			if (escolha == encontrados)
			{
				correto = true;
				break;
			}
		}

		//getchar();

		return correto;
	}

	static vector<vector<int> >& verifica_sol(CUTINFOptr cutinfo, double* sol)
	{
		RepPRVP rep(cutinfo->P);
		vector<vector<int> > posicoes(cutinfo->P);

		for (int i = 0; i < cutinfo->numeroDeRotas; i++)
		{
			if (sol[i] >= 1 - EpsForIntegrality)//rota está na solucao

			{
				posicoes.at(cutinfo->periodos.at(i)).push_back(i);
				rep.at(cutinfo->periodos.at(i)).push_back(cutinfo->rotas.at(i));
			}
		}

		vector<vector<int> >& incompativeis = *new vector<vector<int> > ;

		for (int c = 1; c <= cutinfo->N; c++)
		{
			vector<int> encontrados;
			vector<int> maybe;
			for (int p = (rep.size() - 1); p >= 0; p--)
			{
				for (int r = 0; r < rep.at(p).size(); r++)
				{
					for (int pos = 0; pos < rep.at(p).at(r).size(); pos++)
					{
						if (rep.at(p).at(r).at(pos) == c)
						{
							maybe.push_back(posicoes.at(p).at(r));
							encontrados.push_back(p);
							break;
						}
					}
				}
			}

			if (!SetPart::estaPadrao(cutinfo, encontrados, c))
			{
				incompativeis.push_back(maybe);
			}

		}

		/*if (incompativeis.size() > 0)
		 {
		 cout << "incomp: " << incompativeis << endl << endl;
		 cout << "posicoes: " << posicoes << endl;
		 cout << "rep: " << rep << endl;
		 exit(1);
		 }*/
		return incompativeis;
	}

	static RepPRVP& carregaSolucaoIncumbente(CPXCENVptr env, CPXLPptr model, INCUMBENTINFOptr incumbentinfo, double *x)
	{

		double *solY = new double[incumbentinfo->rotas.size()];
		RepPRVP& rep = *new RepPRVP(incumbentinfo->P);
		char name[100];
		for (int k = 0; k < incumbentinfo->numcols; k++)
		{
			string s = getCollumName(env, model, k);
			int i;
			if (s.at(0) == 'Y')
			{
				strncpy(name, s.c_str(), 100);
				sscanf(name, "Y(%d)", &i);
				solY[i] = x[k];
				if (solY[i] >= 1 - EpsForIntegrality)
				{
					rep.at(incumbentinfo->periodos.at(i)).push_back(incumbentinfo->rotas.at(i));
				}
			}
		}

		delete[] solY;

		for (int p = 0; p < rep.size(); p++)
			if (rep.at(p).size() < incumbentinfo->pPRVP->nVeic)
				rep.at(p).push_back(vector<int> ());

		return rep;
	}

	static int CPXPUBLIC
	myincumbentcheck(CPXCENVptr env, void *cbdata, int wherefrom, void *cbhandle, double objval, double *x, int *isfeas_p, int *useraction_p)
	{

		int iterations;
		CPXgetcallbackinfo(env, cbdata, wherefrom, CPX_CALLBACK_INFO_MIP_ITERATIONS, &iterations);
		int status = 0;

		if (iterations > 0)
		{

			pthread_mutex_lock(&cs_mutex);

			*useraction_p = CPX_CALLBACK_DEFAULT;

			// Instance data


			INCUMBENTINFOptr incumbentinfo = (INCUMBENTINFOptr) cbhandle;

			CPXLPptr model = incumbentinfo->model;
			int numcols = incumbentinfo->numcols;
			vector<vector<int> > rotas = incumbentinfo->rotas;
			vector<int> periodos = incumbentinfo->periodos;
			BasicIteratedLocalSearchPRVP* ils = incumbentinfo->ils;
			vector<vector<int> >& rotasDoCliente = incumbentinfo->rotasDoCliente;

			//int rSize = 0;
			double eps = 0.00001;

			// Vetores e variaveis do ILS

			int melhoras = 1;
			double custoILS = objval;
			int *rotaSizeILS;
			double *custoRotaILS;

			vector<vector<int> > solucaoILS;
			vector<vector<double> > loadILS;

			// Carrega a solucao incumbente

			RepPRVP& rep = SetPart::carregaSolucaoIncumbente(env, model, incumbentinfo, x);

			// Executa a heurística
			PRVPEvaluator ev(*incumbentinfo->pPRVP);
			LPool l(*incumbentinfo->pPRVP, incumbentinfo->tolerancia);
			EvaluationPRVP& e = ev.evaluate(rep);
			SolutionPRVP& sol = *new TestSolution<RepPRVP> (rep);
			ils->exec(sol, e, 0, 0, l);

			custoILS = e.evaluation();

			// Atualiza a melhor solucao

			if (custoILS < custoBestOfAll - eps)
			{
				//atualizando a solução
				solution = &sol;
				evaluation = &e;
				custoBestOfAll = custoILS;
				//Imprimindo resultado
				cout << "Melhora ILS-CPLEX: " << custoBestOfAll << endl;
				//ImprimeResultadoParcial(objval, custoBestOfAll, 0, melhoras, VizILS, Perturbacao);
				cout << endl;
				cutOff = custoBestOfAll;

				double dualval;
				CPXgetcallbackinfo(env, cbdata, wherefrom, CPX_CALLBACK_INFO_BEST_REMAINING, &dualval);
				if (custoBestOfAll < dualval)
				{
					*useraction_p = CPX_CALLBACK_FAIL;
					cout << "\nOtimização abortada! Solucao da heuristica melhor que o Limite Inferior do CPLEX." << endl;
				}
			}

			//getchar();
			delete[] rotaSizeILS;
			delete[] custoRotaILS;
			//acabei de adicionar
			delete &sol;
			delete &e;
			delete &rep;

			pthread_mutex_unlock(&cs_mutex);
		}

		return (status);

	}

	static int CPXPUBLIC
	mycutcallback(CPXCENVptr env, void *cbdata, int wherefrom, void *cbhandle, int *useraction_p)
	{
		int status = 0;

		// Mutex

		pthread_mutex_lock(&cs_mutex);

		*useraction_p = CPX_CALLBACK_DEFAULT;

		// Getting the node's depth

		int depth;
		CPXgetcallbacknodeinfo(env, cbdata, wherefrom, 0, CPX_CALLBACK_INFO_NODE_DEPTH, &depth);

		double dualval;
		CPXgetcallbackinfo(env, cbdata, wherefrom, CPX_CALLBACK_INFO_BEST_REMAINING, &dualval);
		if (custoBestOfAll < dualval)
		{
			*useraction_p = CPX_CALLBACK_FAIL;
			cout << "\nOtimização abortada! Solucao da heuristica melhor que o Limite Inferior do CPLEX." << endl;
		}

		//if (depth <= 7 && *useraction_p == CPX_CALLBACK_DEFAULT)

		{

			//int status = 0;

			// Instance data

			CUTINFOptr cutinfo = (CUTINFOptr) cbhandle;

			CPXLPptr model = cutinfo->model;
			int numcols = cutinfo->numcols;
			int N = cutinfo->N;
			double *x = cutinfo->x;
			double Q = cutinfo->Q;
			double *delivery = cutinfo->delivery;
			int numeroDeRotas = cutinfo->numeroDeRotas;
			vector<vector<int> > rotasDoCliente = cutinfo->rotasDoCliente;
			vector<vector<int> > rotas = cutinfo->rotas;

			// Vehicle Capacity (ceil)

			int CAP = ceil(Q);

			// Number of customers

			int NoOfCustomers = N;

			// Sum of the delivery and pickup loads

			int sumDelivery = 0;

			// Demands (Demand -> delivery; Demand2 -> pickup)

			int *Demand = new int[N + 1];

			for (int i = 1; i < N + 1; i++)
			{
				Demand[i] = floor(delivery[i]);
				sumDelivery += Demand[i];
			}

			// Allocationg memory for the cuts

			int CustListSize = 0;
			int *CustList = new int[N + 1];

			// Array that stores the fractional solution

			double *sol = new double[numcols + 1];
			for (int i = 0; i < numcols + 1; i++)
			{
				sol[i] = 0;
			}

			// Loading the fractional solution and storing in the array sol

			CPXgetcallbacknodex(env, cbdata, wherefrom, x, 0, numcols - 1);
			SetPart::storeLPSolution(env, model, numcols, x, sol);

			/*for (int i = 0; i < cutinfo->numeroDeRotas; i++)
			 {
			 if (sol[i] != 0)//rota está na solucao
			 {
			 cout<<"sol: "<<sol[i]<<endl;
			 cout<<"["<<cutinfo->periodos.at(i)<<"]"<<cutinfo->rotas.at(i)<<endl;
			 }
			 }*/

			// Computing the number of edges (non-zeros)

			// Flag that check if there is at least one variable with fractional cost
			bool isFractional = SetPart::isFrac(numcols, sol);

			if (!isFractional)
			{
				//cout << "solucao inteira" << endl;
				//getchar();
				// VALIDAÇÃO DA SOLUÇÃO PARA O PROBLEMA
				vector<vector<int> >& rotas_incompativeis = SetPart::verifica_sol(cutinfo, sol);
				// guardar o número das rotas que forem inviáveis : nome vector: rotas_incompativeis
				if (rotas_incompativeis.size() > 0)
				{
					int colindex;
					char objName[100];

					for (int r = 0; r < rotas_incompativeis.size(); r++)
					{
						vector<int> cutindAux;
						for (int i = 0; i < rotas_incompativeis.at(r).size(); i++)
						{
							sprintf(objName, "Y(%d)", rotas_incompativeis.at(r).at(i));
							CPXgetcolindex(env, model, objName, &colindex);
							cutindAux.push_back(colindex);
						}

						int cutnz = cutindAux.size();
						int *cutind = new int[cutnz];
						double *cutcoef = new double[cutnz];

						for (int i = 0; i < cutnz; i++)
						{
							cutind[i] = cutindAux[i];
							cutcoef[i] = 1;
						}
						double RHS = rotas_incompativeis.at(r).size() - 1;
						status = CPXcutcallbackadd(env, cbdata, wherefrom, cutnz, RHS, 'L', cutind, cutcoef, 1);
						if (status)
						{
							fprintf(stderr, "Failed to add cut.\n");
						}
						else
						{
							cout << "PRVP Cut applied" << endl;
							//numCutsMB++;
						}
						delete[] cutind;
						delete[] cutcoef;
					}
					//acabei de adicionar
					delete &rotas_incompativeis;

				}
				else
				{

				}

			}

			delete[] CustList;
			delete[] Demand;

			delete[] sol;
		}

		pthread_mutex_unlock(&cs_mutex);
		//cout << "###rodou!!!" << endl;
		return status;
	}

	// Function isFrac

	static bool isFrac(int numcols, double *sol)
	{

		for (int i = 0; i < numcols; i++)
		{
			if ((sol[i] > EpsForIntegrality && sol[i] < 1 - EpsForIntegrality))
			{
				return true;
			}
		}

		return false;
	}

public:

	SetPart(PRVPProblemInstance& _pPRVP, EvaluationPRVP& _e, SolutionPRVP& _sol) :
		pPRVP(_pPRVP)
	{
		custoBestOfAll = _e.evaluation();
		evaluation = &_e;
		solution = &_sol;
		cutOff = 99999999;//infinito
	}
	;

	~SetPart()
	{
	}

	void geraMST(LPool& lista, SolutionPRVP& s)
	{
		int num_rotas = lista.rotas.size();

		int num_rotas_validas = 0;

		RepPRVP& rep = s.getR();

		for (int p = 0; p < rep.size(); p++)
		{
			for (int r = 0; r < rep.at(p).size(); r++)
			{
				if (rep.at(p).at(r).size() != 0)
					num_rotas_validas++;
			}
		}

		int totalDeRotas = num_rotas + num_rotas_validas;

		int *mipStart = new int[totalDeRotas];
		for (int i = 0; i < num_rotas; i++)
		{
			mipStart[i] = 0;
		}
		for (int i = num_rotas; i < num_rotas + num_rotas_validas; i++)
		{
			mipStart[i] = 1;
		}

		ofstream mst("prvp.mst", ios::out);
		if (!mst)
		{
			cerr << "Arquivo nao pode ser aberto" << endl;
		}
		mst << "NAME" << endl;
		for (int i = 0; i < num_rotas + num_rotas_validas; i++)
		{
			mst << "    Y(" << i << ") " << mipStart[i] << endl;
		}
		mst << "ENDATA";
		mst.close();

		//adicionando a solução corrente a lista
		lista.add_withoutTest(s);

		delete[] mipStart;
	}

	void geraLP(LPool& lista)
	{
		int numeroDeRotas = lista.rotas.size();
		int N = pPRVP.nCli;

		// Inicia o ambiente Concert
		IloEnv env;
		IloModel modelo(env);

		// Variaveis y (correspondem as rotas)

		IloBoolVarArray y(env, numeroDeRotas);
		char var[100];
		for (int i = 0; i < numeroDeRotas; i++)
		{
			sprintf(var, "Y(%d)", (int) i);
			y[i].setName(var);
			modelo.add(y[i]);
		}

		// Funcao Objetivo

		IloExpr FO(env);
		for (int i = 0; i < numeroDeRotas; i++)
		{
			FO += lista.custo.at(i) * y[i];
		}

		modelo.add(IloMinimize(env, FO));

		// Restricoes se particionamento de conjunto

		for (int i = 1; i < N + 1; i++)
		{
			IloExpr temp(env);
			for (int j = 0; j < lista.cliInRoute->at(i).size(); j++)
			{
				temp += y[lista.cliInRoute->at(i).at(j)];
			}
			modelo.add(temp == pPRVP.dados->at(i).at(NVISITAS));
		}

		// Restricoes de cliente visitado no máximo uma vez por periodo

		for (int i = 1; i < N + 1; i++)
		{
			for (int p = 0; p < pPRVP.nPeriod; p++)
			{
				IloExpr temp(env);
				for (int j = 0; j < lista.cliInRoute->at(i).size(); j++)
				{
					if (p == lista.periodo.at(lista.cliInRoute->at(i).at(j)))
					{
						temp += y[lista.cliInRoute->at(i).at(j)];
					}
				}
				modelo.add(temp <= 1);
			}
		}

		// Restricoes que indicam o numero de rotas da solucao


		for (int p = 0; p < pPRVP.nPeriod; p++)
		{
			IloExpr temp(env);
			for (int r = 0; r < lista.periodo.size(); r++)
			{
				if (p == lista.periodo.at(r))
				{
					temp += y[r];
				}
			}
			modelo.add(temp <= pPRVP.nVeic);
		}

		// Gera LP

		IloCplex PRVP(modelo);
		PRVP.exportModel("lp.lp");

		//Resolver direto pelo Concert
		//PRVP.setParam(IloCplex::PreLinear, 0);
		//PRVP.readMIPStart("prvp.mst");
		//PRVP.setParam(IloCplex::AdvInd, 1); //AdvInd = 1 ou 2

		//PRVP.solve();

		PRVP.end();
		env.end();

	}

	void ambienteCPlex(LPool& lista, BasicIteratedLocalSearchPRVP& ils)
	{
		Timer t;
		// Variaveis auxiliares do ambiente cplex
		double before = t.now();
		int status = 0;
		int solstat, cur_numcols;

		// Inicializa o ambiente cplex

		CPXENVptr env = NULL;
		CPXLPptr model = NULL;

		env = CPXopenCPLEX(&status);
		model = CPXcreateprob(env, &status, "lp");
		status = CPXreadcopyprob(env, model, "lp.lp", NULL);
		cur_numcols = CPXgetnumcols(env, model);

		// Parametros do cplex

		CPXsetintparam(env, CPX_PARAM_MIPINTERVAL, 100);
		CPXsetintparam(env, CPX_PARAM_SCRIND, CPX_ON);
		CPXsetintparam(env, CPX_PARAM_PRELINEAR, 0);
		CPXsetintparam(env, CPX_PARAM_MIPCBREDLP, CPX_OFF);
		CPXsetintparam(env, CPX_PARAM_THREADS, 1);
		CPXsetintparam(env, CPX_PARAM_REDUCE, CPX_PREREDUCE_PRIMALONLY);
		//CPXsetintparam(env, CPX_PARAM_MIPEMPHASIS, 0);
		//CPXsetintparam(env, CPX_PARAM_VARSEL, 3);

		// Informacoes necessarias para a busca local na solucao incumbente
		//cout << "antes" << endl;
		INCUMBENTINFO incumbentinfo;
		atualizaIncumbentinfo(&incumbentinfo, model, ils, lista);

		// Informacoes necessarias para geracao de cortes

		CUTINFO cutinfo;
		atualizaCutinfo(&cutinfo, model, lista);

		// Ativando os callbacks

		//===CPXsetbranchcallbackfunc(env, usersetbranch, NULL);
		CPXsetcutcallbackfunc(env, mycutcallback, &cutinfo);
		CPXsetincumbentcallbackfunc(env, myincumbentcheck, &incumbentinfo);

		// Carrega solucao inicial

		CPXreadcopymipstarts(env, model, "prvp.mst");

		// Imprime tempo da heuristica + pre-processamento da fase hibrida

		double meio = t.now();
		cout << "\nTempo da heuristica + pre-processamento da fase hibrida: " << meio - before << "\n" << endl;

		double cplexTimeBefore, cplexTimeAfter;

		// Otimiza

		CPXgettime(env, &cplexTimeBefore);
		CPXmipopt(env, model);
		CPXgettime(env, &cplexTimeAfter);

		// Imprime o tempo de cpu do cplex

		cout << "\nTempo do exato: " << cplexTimeAfter - cplexTimeBefore << endl;

		double after = (meio - before) + (cplexTimeAfter - cplexTimeBefore);

		// Atualiza a solucao final

		//===AtualizaSolucaoFinal(env, model, cur_numcols, numeroDeRotas, &custoBestOfAll, &rSizeBestOfAll, rotaSizeBestOfAll, solucaoBestOfAll, rotas, rSize, numVeicDepositos, depotDaRota, eps);
		//===rSizeBestOfAll = rSize;

		// Fecha o ambiente cplex

		CPXfreeprob(env, &model);
		CPXcloseCPLEX(&env);
	}

	SolutionPRVP& getSolution()
	{
		return *solution;
	}

	EvaluationPRVP& getEvaluation()
	{
		return *evaluation;
	}

};

#endif /*SETPART_HPP*/

